/*
Previous Scheme

   1. An initial clustering is performed using full length alignments (60% of longer)
   2. All clusters with > 10 members are saved as seed clusters
   3. Then fragment (non-full-length) alignments are included (60% of longer OR 75% of 
		shorter)
   4. Members are parachuted into the initial clusters in a non-transitive way.
      I.e. only proteins with direct links to the seed clusters are included.
   5. Then the remaining unclustered proteins are clustered. This explain why seed 
		clusters need to be > 10. So that small clusters can be joined together using 
		fragments.

Bug ALERT! ALERT!

   1. Because of the privileged state of seed clusters, any related protein will be 
			parachuted into a seed cluster and unavailable for clustering otherwise.
   2. This is especially serious because
         1. Full length alignment in small clusters are excluded from seed clusters
         2. Some of the proteins in full length alignments will have poor partial alignmen
			ts to the seed clusters

Solution

   1. Parachute transitively and recursively
   2. Do not exclude parachuted proteins from further clustering

	Solution 2 was chosen because we have experience with non-transitive 
	clustering and because this would require the least amount of code 
	modification. (Andreas uses solution 1) 

*/
#include <string>
#include <sstream>
#include <istream>
#include <algorithm>
#include <numeric>
#include <istream>
#include <limits>
#include <vector>
#include <fstream>
#include <deque>
#include <map>
#include <stdexcept>
#include <tr1/unordered_map>


// boost functions
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/foreach.hpp>


// my functions
#include <stlhelper.h>
#include <tab_delimited_str.h>
#include <single_linkage.h>
#include <portable_timer.h>
#include <progress_indicator.h>
#include <print_error.h>
#include <ultoa.h>
#include <index_string.h>
#define DEBUG_LEVEL 5
#include <seq_pairs.h>
#include <debug_func.h>

#include "get_arg.h"

using std::lower_bound;
using std::partial_sort_copy;
using std::runtime_error;
using std::cerr;
using std::istream;
using std::ifstream;
using std::ofstream;
using std::istringstream;
using std::cout;
using std::cin;
using std::string;
using std::vector;
using std::deque;
using std::sort;
using boost::lexical_cast;
using boost::algorithm::replace_all;
using std::ostringstream;
using stlhelper::lookup;





typedef deque<vector<unsigned> >	deque_vec_u;
typedef std::tr1::unordered_map<unsigned int, unsigned int> t_hashuu;

struct t_clustered_id_pair : public std::pair<unsigned, unsigned>
{
	float dist;
	unsigned cluster_id;
	t_clustered_id_pair(unsigned id1_, unsigned id2_, float dist_)
		:std::pair<unsigned, unsigned>(id1_, id2_), dist(dist_), cluster_id(0)
	{
	}
};

typedef deque<t_clustered_id_pair> t_pairs;

bool operator < (const t_clustered_id_pair& a, const t_clustered_id_pair& b)
{
	return a.cluster_id < b.cluster_id;
}

//________________________________________________________________________________________

//		parse_blast_pairs_from_file

//________________________________________________________________________________________
struct do_parse_pairs
{
	t_index_string&			prot_id_str;
	t_pairs& 				prot_pairs;
	std::string				comment_letters;
	do_parse_pairs(t_index_string&	prot_ids_, t_pairs& prot_pairs_, 
				   const std::string&	comment_letters_):
				prot_id_str(prot_ids_),  prot_pairs(prot_pairs_), 
				comment_letters(comment_letters_){}

	void operator()(istream& is, unsigned)
	{
		// get data for one protein pair
		string name1, name2;

		// comment if false
		string comment_line;
		if (!seq_pairs::parse_ids(is, name1, name2, comment_line, comment_letters,
								  seq_pairs::ePARSE_IDS_ONLY))
			return;


		float dist;
		is >> dist;
		if (!is)
			throw runtime_error("The pairs distances data is not in the correct format");

		// ignore any other fields
		is.ignore(32768, '\n');
		is.putback('\n');
		

		// get protein and gene indices
		prot_pairs.push_back(t_clustered_id_pair(prot_id_str.index(name1),  
												 prot_id_str.index(name2), dist));
	}
};


//________________________________________________________________________________________

//		parse_blast_pairs_from_file

// 				offset is the index of the first cluster
//				= num of predefined clusters
//________________________________________________________________________________________
void create_prot_ids_to_cluster_lookup(	const deque_vec_u& 	prot_id_clusters,
										vector<unsigned>& 	prot_id_to_cluster_id,
										unsigned			offset)
{
	debug_func_cerr(6);
	// cluster_id 0 reserved for orphans
	// increment no longer necessary because cluster_id now points to a real cluster #0
	//++offset; 
	
	for (unsigned ii = 0; ii < prot_id_clusters.size(); ++ii)
		for (unsigned jj = 0; jj < prot_id_clusters[ii].size(); ++jj)
			prot_id_to_cluster_id[prot_id_clusters[ii][jj]] = ii + offset;
}




//________________________________________________________________________________________

//		assign_unclustered_prot_pairs_to_clusters

// 				Lookup cluster of protein pairs which have not yet been
//				clustered
//________________________________________________________________________________________
void assign_unclustered_prot_pairs_to_clusters(	t_pairs&					prot_pairs,
												const vector<unsigned>& 	prot_id_to_cluster_id)
{
	debug_func_cerr(6);
	for (unsigned ii = 0; ii < prot_pairs.size(); ++ii)
		if (!prot_pairs[ii].cluster_id)
			prot_pairs[ii].cluster_id = prot_id_to_cluster_id[prot_pairs[ii].first];
}


//________________________________________________________________________________________

//		mark_prot_pairs_of_small_clusters_for_removal

//________________________________________________________________________________________
void mark_prot_pairs_of_small_clusters_for_removal(	t_pairs&				prot_pairs,
													const vector<unsigned>&	cluster_id_to_cluster_size)
{
	debug_func_cerr(6);
	unsigned to_be_removed = std::numeric_limits<unsigned>::max();
	for (unsigned ii = 0; ii < prot_pairs.size(); ++ii)
		if (!cluster_id_to_cluster_size[prot_pairs[ii].cluster_id])
			prot_pairs[ii].cluster_id = to_be_removed;
}


//________________________________________________________________________________________

//		index_of_first_clustered_pair

// 				index of first pair whose cluster_id is > 0
// 				pairs should already be sorted order (by cluster_id)
// 
//________________________________________________________________________________________
unsigned index_of_first_clustered_pair(const t_pairs& prot_pairs)
{
	debug_func_cerr(6);
	t_clustered_id_pair first_clustered_pair(0,0,0.0);
	first_clustered_pair.cluster_id = 1;
	return lower_bound(prot_pairs.begin(), prot_pairs.end(), 
								first_clustered_pair) - prot_pairs.begin();
}

//________________________________________________________________________________________

//		parachute_prot_into_cluster

//________________________________________________________________________________________
void parachute_prot_into_cluster(unsigned 		prot_id, 
								 unsigned		cluster_id,
								 vector<char>&	is_parachuted,
								 t_hashuu&		prot_id_to_parachute_id,
								 deque_vec_u&	parachute_id_to_cluster_ids,
								 deque_vec_u&	prot_id_clusters)
{
	debug_func_cerr(6);
	// if prot_id has not previously parachuted into any cluster: 
	// 	created parachuted data structures and do it now
	if (!is_parachuted[prot_id])
	{
		is_parachuted[prot_id] = 1;
		prot_id_to_parachute_id[prot_id] = parachute_id_to_cluster_ids.size();
		parachute_id_to_cluster_ids.push_back(vector<unsigned>());
		parachute_id_to_cluster_ids.back().push_back(cluster_id);
		prot_id_clusters[cluster_id].push_back(prot_id);
		debug_cerr(6, "parachuted " << prot_id << " into " << cluster_id);
		return;
	}

	// if already has been parachuted,
	// 		was it this cluster which it has been parachuted into?
	
	unsigned parachute_id = prot_id_to_parachute_id[prot_id];
	if (!includes(parachute_id_to_cluster_ids[parachute_id], cluster_id))
	{
		prot_id_clusters[cluster_id].push_back(prot_id);
		parachute_id_to_cluster_ids[parachute_id].push_back(cluster_id);
		debug_cerr(6, "parachuted " << prot_id << " also into " << cluster_id);
	}
}



//________________________________________________________________________________________

//		get_common_items

//			get vector items in common

//			assumes both vectors are very small:
// 			therefore don't need to presort

//________________________________________________________________________________________
template<class CONT >
void get_common_items(CONT& in_common, 
						const CONT& array1,
						const CONT& array2)
{
	debug_func_cerr(6);
	typedef typename CONT::const_iterator citer;
	citer i = array1.begin();
	citer end = array1.end();
	while (i != end)
	{
		if (array2.end() != find(array2.begin(), array2.end(), *i) )
			in_common.push_back(*i);
		++i;
	}
}

//________________________________________________________________________________________

//		index_name_of_multiply_parachuted
// 
//		 parachuted proteins may be in multiple clusters
//		 	add index to name to indicate this 
 

//________________________________________________________________________________________
void index_name_of_multiply_parachuted(	string& 			original_prot_id_str,
										const vector<char>&	is_parachuted,
										const t_hashuu&		prot_id_to_parachute_id,
										const deque_vec_u&	parachute_id_to_cluster_ids,
										unsigned			prot_id,
										unsigned			cluster_id,
										const string&		substitution_format)
{
	debug_func_cerr(6);
	unsigned parachute_id = lookup(prot_id).in(prot_id_to_parachute_id);
	if (parachute_id_to_cluster_ids[parachute_id].size() > 1)
	{
		// The cluster_id may not be found in parachute_id_to_cluster_ids if
		//		this cluster has been formed by single linkage AFTER and separately
		//		from parachuting.
		// In this case, the index will be the number of parachuted clusters (which is OK)
		string index = ultoa(find_pos_in_vector(parachute_id_to_cluster_ids[parachute_id], 
												cluster_id));
		string new_prot_id_str = substitution_format;
		replace_all(new_prot_id_str, "[ID]", original_prot_id_str);
		replace_all(new_prot_id_str, "[INDEX]", index);
		original_prot_id_str = new_prot_id_str;
	}										
}



//________________________________________________________________________________________

//		resolve_parachuted_pairs
// 
// 			Remove all pairs which have a mix of parachuted & non-parachuted prot_ids
// 						i.e. parachuting is not a transitive process!
// 			Keep pairs which have been parachuted into the same clusters
// 						If pairs share > 1 cluster, 
// 						each cluster should have its own copy of the pair			
// 
//________________________________________________________________________________________
void resolve_parachuted_pairs(	t_pairs& 			prot_pairs,
								const deque_vec_u&	parachute_id_to_cluster_ids,
								const vector<char>&	is_parachuted,
								const t_hashuu&		prot_id_to_parachute_id,
								t_program_args&		args)
{
	args.VERBOSE(2) << "\tRemoving pairs with mixture of parachuted and unclustered...\n";
	args.VERBOSE(3) << "\t\t(Parachuting is not a transitive process)\n";

	unsigned cnt_pairs_one_parachuted = 0;
	unsigned cnt_pairs_parachuted_different_clusters = 0;

	// for pairs parachuted into multiple clusters
	t_pairs				extra_prot_pairs;
	// mark mixed pairs
	BOOST_FOREACH(t_clustered_id_pair& id_pair, prot_pairs) 
	{
		if (id_pair.cluster_id)
			continue;

		// if only one out of pair parachuted, to be single linkage clustered later
		if (is_parachuted[id_pair.first] !=
			is_parachuted[id_pair.second])
		{	
			cnt_pairs_one_parachuted++;
			continue;
		}

		// neither parachuted: to be single linkage clustered later
		if (!is_parachuted[id_pair.first])
			continue;

		// 
		// both parachuted:
		// 		Any clusters in common?
		// 		Each cluster should have its own copy of the pair
		// 		Add extras in as duplicate
		// 
		vector<unsigned> common_cluster_ids;
		unsigned parachute_id1 = lookup(id_pair.first ).in(prot_id_to_parachute_id);
		unsigned parachute_id2 = lookup(id_pair.second).in(prot_id_to_parachute_id);
		get_common_items(	common_cluster_ids, 
							parachute_id_to_cluster_ids[parachute_id1],
							parachute_id_to_cluster_ids[parachute_id2]);

		// No: parachuted into different clusters
		if (common_cluster_ids.empty())
		{	
			++cnt_pairs_parachuted_different_clusters;
			continue;
		}

		// each cluster in common should have its own copy of the pair
		BOOST_FOREACH(unsigned i, common_cluster_ids) 
		{
			extra_prot_pairs.push_back(id_pair);
			extra_prot_pairs.back().cluster_id = i;
		}
	}
	

	args.VERBOSE(3)	<< "\t\t" << cnt_pairs_one_parachuted
					<< " pairs clustered later because only one has been parachuted\n";
	args.VERBOSE(3)	<< "\t\t" << cnt_pairs_parachuted_different_clusters
					<< " pairs clustered later because they have been "
						"parachuted into different clusters\n";
	args.VERBOSE(3)	<< "\t\t" << extra_prot_pairs.size()
					<< " pairs added because only they are present in multiple clusters\n";

	// add protein pairs in multiply copies
	prot_pairs.insert(prot_pairs.end(), extra_prot_pairs.begin(),  extra_prot_pairs.end());
}

//________________________________________________________________________________________

//		removed_marked_pairs

//			remove pairs marked for deletion
//			N.B. Remember to sort first

//________________________________________________________________________________________
void removed_marked_pairs(t_pairs&				prot_pairs,
						  t_program_args&		args)
{
	
	// create pair which matches the "deletion" cluster
	t_clustered_id_pair pair_to_remove(0,0,0.0);
	unsigned mark_for_removal = std::numeric_limits<unsigned>::max();
	pair_to_remove.cluster_id = mark_for_removal;

	
	t_pairs::iterator removal_beg = lower_bound(	prot_pairs.begin(), 
													prot_pairs.end(), 
													pair_to_remove);
	if (removal_beg != prot_pairs.end())
	{
		args.VERBOSE(3) << "\t" << (prot_pairs.end() - removal_beg)
									<< " pairs removed in total.\n";
		prot_pairs.erase(removal_beg, prot_pairs.end());
	}
	
}
//________________________________________________________________________________________

//		cluster_by_predefined_clusters
// 
// 			A) Both pairs belong to the same predefined cluster :
// 						cluster them there
// 			B) Both pairs belong to different predefined clusters :
// 						mark for deletion.
// 						No cross joining
// 			C) Only one belongs to a predefined cluster :
// 						Parachute it
// 
//________________________________________________________________________________________
void cluster_by_predefined_clusters(t_pairs&				prot_pairs,
									const t_index_string&	prot_id_strs,
									const vector<unsigned>&	prot_id_to_cluster_id,
									vector<char>&			is_parachuted,
									t_hashuu&				prot_id_to_parachute_id,
									deque_vec_u&			parachute_id_to_cluster_ids,
									deque_vec_u& 			prot_id_clusters,
									t_program_args& 		args)
{
	args.VERBOSE(2) << "\tCluster protein pairs by predefined clusters...\n";


	// This cluster_id is for unwanted pairs where both are already in different
	// 	predefined clusters.
	unsigned mark_for_removal = std::numeric_limits<unsigned>::max();

	
	// pairs cross linking different clusters are discarded
	unsigned cnt_pairs_different_clusters = 0;

	BOOST_FOREACH(t_clustered_id_pair& id_pair, prot_pairs) 
	{
		unsigned cluster_id1 = prot_id_to_cluster_id[id_pair.first];
		unsigned cluster_id2 = prot_id_to_cluster_id[id_pair.second];


		//  both zero 		= unclustered
		//  both non-zero	= in same cluster
		if (cluster_id1 == cluster_id2)
		{
			id_pair.cluster_id = cluster_id1;
			continue;
		}

		//	both non-zero & different = cross-join clustered. Delete
		if (cluster_id1 && cluster_id2)
		{
			id_pair.cluster_id = mark_for_removal;
			++cnt_pairs_different_clusters;
			continue;
		}

		//	one zero: parachute unclustered (zero) into clustered
		//	id_pair.cluster_id is the cluster_id + 1 (cluster = 0 when unclustered)
		if (cluster_id1)
		{
			id_pair.cluster_id = cluster_id1;

			parachute_prot_into_cluster(id_pair.second, 
										id_pair.cluster_id,
										is_parachuted,
										prot_id_to_parachute_id,
										parachute_id_to_cluster_ids,
										prot_id_clusters);
		}
		else
		{
			id_pair.cluster_id = cluster_id2;

			parachute_prot_into_cluster(id_pair.first, 
										id_pair.cluster_id,
										is_parachuted,
										prot_id_to_parachute_id,
										parachute_id_to_cluster_ids,
										prot_id_clusters);
		}
	}


		

	args.VERBOSE(3) << "\t\t" << prot_id_to_parachute_id.size()
						<< " identifiers parachuted into pre-existing clusters\n";
	args.VERBOSE(3) << "\t\t" << cnt_pairs_different_clusters
						<< " pairs removed because both were already in "
						   "distinct pre-existing clusters\n";
}

//________________________________________________________________________________________

//		dump_id_pairs

//			for debugging

//________________________________________________________________________________________
#ifndef NDEBUG
void dump_id_pairs(	const t_pairs& 			prot_pairs_orig, 
					const t_index_string&	prot_id_strs)
{
	t_pairs prot_pairs = prot_pairs_orig;
	sort(prot_pairs.begin(), prot_pairs.end());
	unsigned curr_cluster = 500000;
	BOOST_FOREACH(t_clustered_id_pair& prot_pair, prot_pairs) 
	{
		if (prot_pair.cluster_id != curr_cluster)
		{
			curr_cluster = prot_pair.cluster_id;
			cerr << ">cluster# " << curr_cluster << "\n";
		}
		cerr	<< prot_id_strs[prot_pair.first]	<< "\t"
				<< prot_id_strs[prot_pair.second]	<< "\t"
				<< prot_pair.dist				<< "\n";
	}
}
#endif


//________________________________________________________________________________________

// Single linkage cluster all pairs not in predefined clusters

//________________________________________________________________________________________
void single_linkage_cluster(const t_pairs&		prot_pairs,
							deque_vec_u&		prot_id_clusters,
							vector<unsigned>&	prot_id_to_cluster_id,
							t_program_args&		args)
{	
	unsigned cnt_unclustered_pairs = index_of_first_clustered_pair(prot_pairs);
	if (cnt_unclustered_pairs == 0)
		return;
	// find (one past) last item to be clustered
	t_pairs::const_iterator unclustered_pairs_end = 
							prot_pairs.begin() + cnt_unclustered_pairs;

	args.VERBOSE(2) << "\tSingle linkage clustering...\n";
	unsigned DOT_SIZE = (cnt_unclustered_pairs / 50) + 1;
	args.VERBOSE(3) << "\t\t(" << DOT_SIZE * 50 << " pairs per line...)\n";
	t_progress_indicator dots(args.VERBOSE(2), DOT_SIZE);


	deque_vec_u	new_prot_id_clusters;
	t_pairs::const_iterator beg = prot_pairs.begin();
	single_linkage_clusters(beg,  unclustered_pairs_end, new_prot_id_clusters, dots);
	dots.finish();

	//
	//	Remember which prot id is which cluster
	//
	create_prot_ids_to_cluster_lookup(new_prot_id_clusters, prot_id_to_cluster_id, prot_id_clusters.size());
	args.VERBOSE(2) << "\t" << new_prot_id_clusters.size() << " clusters "
														"created by single linkage...\n";


	// add to predefined clusters
	for (unsigned i = 0; i < new_prot_id_clusters.size(); ++i)
	{
		prot_id_clusters.push_back(std::vector<unsigned>());
		prot_id_clusters.back().swap(new_prot_id_clusters[i]);
	}
}

//________________________________________________________________________________________

//		print_identifiers

//________________________________________________________________________________________
template <typename CONT>
struct op_cmp_by_lookup
{
	const CONT& container;
	op_cmp_by_lookup(const CONT& cont) :container(cont){}
	bool operator()(unsigned i,  unsigned j) const
	{
		return container[i] < container[j];
	}
};
void sort_and_print_identifiers(	const t_index_string&	prot_id_strs,
									deque_vec_u&			prot_id_clusters,
									const vector<unsigned>&	cluster_id_to_cluster_size,
									t_program_args&			args)

{
	args.VERBOSE(2) << "\tOutput cluster of identifiers...\n";
	t_progress_indicator dots(args.VERBOSE(2),  100);
	op_cmp_by_lookup<t_index_string> sort_by_prot_id(prot_id_strs);
	// cluster 0 is for unclustered
	for (unsigned i = 1; i < prot_id_clusters.size(); ++i)
	{
		++dots;
		// ignore clusters which are too small
		if (!cluster_id_to_cluster_size[i])
			continue;
		sort(prot_id_clusters[i].begin(), prot_id_clusters[i].end(), sort_by_prot_id);
		cout << prot_id_strs.string_at(prot_id_clusters[i][0]);
		for (unsigned jj = 1; jj < prot_id_clusters[i].size(); ++jj)
			cout << "\t" << prot_id_strs.string_at(prot_id_clusters[i][jj]);
		cout << "\n";
	}
}


unsigned get_identifier_count(	const t_index_string&	prot_id_strs,
								deque_vec_u&			prot_id_clusters,
								const vector<unsigned>&	cluster_id_to_cluster_size,
								t_program_args&			args)
{
	args.VERBOSE(2) << "\tGet count of identifiers...\n";
	t_progress_indicator dots(args.VERBOSE(2),  100);
	std::vector<unsigned> all_above_threshold_prot_ids;
	
	// cluster 0 is for unclustered
	for (unsigned i = 1; i < prot_id_clusters.size(); ++i)
	{
		++dots;
		// ignore clusters which are too small
		if (!cluster_id_to_cluster_size[i])
			continue;
		all_above_threshold_prot_ids.insert(all_above_threshold_prot_ids.end(), prot_id_clusters[i].begin(), prot_id_clusters[i].end());
	}
	std::sort(all_above_threshold_prot_ids.begin(), all_above_threshold_prot_ids.end());
	remove_duplicates(all_above_threshold_prot_ids);
	return all_above_threshold_prot_ids.size();
}



//________________________________________________________________________________________

//		print_pairs

//________________________________________________________________________________________
void print_pairs(   const t_pairs&			prot_pairs,
					const t_index_string&	prot_id_strs,
					const vector<char>&		is_parachuted,
					const t_hashuu&			prot_id_to_parachute_id,
					const deque_vec_u&		parachute_id_to_cluster_ids,
					const vector<unsigned>&	cluster_id_to_cluster_size,
					t_program_args&			args)
{
	// 
	// find first non-orphan
	//
	unsigned cnt_of_orphans = index_of_first_clustered_pair(prot_pairs);
	if (cnt_of_orphans)
		args.VERBOSE(3) << "\t" << cnt_of_orphans << " self-self pairs were discarded!\n";

	args.VERBOSE(2) << "\tOutput clusters...\n";
	unsigned DOT_SIZE = (prot_pairs.size() + 100 - cnt_of_orphans) / 100;
	args.VERBOSE(3) << "\t\t(dots in two lines...)\n";
	t_progress_indicator 				dots(args.VERBOSE(2),  DOT_SIZE);
	dots.use_timer();

	unsigned print_cluster_id = 0;
	//	current index into clusters etc.
	unsigned curr_cluster_id = std::numeric_limits<unsigned>::max();
	bool distinguish_instances = args.distinguish_instances != "[ID]";
	for (unsigned i = cnt_of_orphans; i < prot_pairs.size(); ++i)
	{
		++dots;
		if (prot_pairs[i].cluster_id != curr_cluster_id)
		{
			curr_cluster_id = prot_pairs[i].cluster_id;
			++print_cluster_id;
			cout	<< ">cluster# "	<< print_cluster_id
					<< " size= "	<< cluster_id_to_cluster_size[curr_cluster_id] << "\n";
		}

		unsigned prot_id1 = prot_pairs[i].first;
		unsigned prot_id2 = prot_pairs[i].second;
		string prot_id1_str = prot_id_strs.string_at(prot_id1);
		string prot_id2_str = prot_id_strs.string_at(prot_id2);
		if (distinguish_instances)
		{
			// parachuted proteins may be in multiple clusters
			// 	add index to name to indicate this 
			if (is_parachuted[prot_id1])
			{
				index_name_of_multiply_parachuted(	prot_id1_str,
													is_parachuted,
													prot_id_to_parachute_id,
													parachute_id_to_cluster_ids,
													prot_id1,
													curr_cluster_id,
													args.distinguish_instances);

			}
			if (is_parachuted[prot_id2])
			{
				index_name_of_multiply_parachuted(	prot_id2_str,
													is_parachuted,
													prot_id_to_parachute_id,
													parachute_id_to_cluster_ids,
													prot_id2,
													curr_cluster_id,
													args.distinguish_instances);

			}
		}
		cout	<< prot_id1_str			<< "\t"
				<< prot_id2_str			<< "\t"
				<< prot_pairs[i].dist	<< "\n";

	}
}

//________________________________________________________________________________________

//		print_summary

//________________________________________________________________________________________
void print_summary(const t_index_string&	prot_id_strs,
				   const t_pairs&			prot_pairs,
				   const deque_vec_u&		prot_id_clusters,
				   const vector<unsigned>&	top_cluster_sizes,
				   unsigned 				cnt_output_clusters,
				   const string&			time,
				   unsigned 				cnt_above_threshold_identifiers,
				   ostream&					strm)
{
	strm  << "\tClustering:\n";
	strm  << "\t" << prot_id_strs.size() << " identifiers in \n";
	strm  << "\t" << prot_pairs.size() << " pairs into \n";
	strm  << "\t" << (prot_id_clusters.size() - 1) << " clusters \n";
	if (cnt_output_clusters != prot_id_clusters.size() - 1)
		strm  << "\t" << cnt_output_clusters << " clusters above threshold size \n";
	if (cnt_above_threshold_identifiers != prot_id_strs.size())
		strm  << "\t" << cnt_above_threshold_identifiers << " identifiers clustered\n";
	strm  << "\ttook " << time << "\n";
	strm  << "\tThe top " << top_cluster_sizes.size() << " cluster sizes are as follows:\n";
	BOOST_FOREACH(unsigned i, top_cluster_sizes) 
		strm << "\t\t" << i << "\n";
	strm  << "\tFinished\n";
	strm  << "\t" << string(60, '_') << "\n\n";
}

//________________________________________________________________________________________

//		read_predefined_clusters

//________________________________________________________________________________________
void read_predefined_clusters(	const t_index_string&	prot_id_strs,
								deque_vec_u&			prot_id_clusters,
								vector<unsigned>&		prot_id_to_cluster_id,
								t_program_args& 		args)
{
	debug_func_cerr(6);
	args.VERBOSE(2) << "\tReading predefined clusters...\n";

	
	ifstream ifs(args.file_cluster.c_str());
	if (!ifs)
		throw std::runtime_error("Could not open " + args.file_cluster);
	unsigned cnt_lines = 0;
	while (ifs)
	{
		string line;
		getline(ifs, line);
		++cnt_lines;
	
		// ignore blank lines
		if (!line.length())
			continue;
	
		//
		// save each index in its appropriate cluster and add cluster lookup
		//
		vector<string> fields;
		boost::split(fields, line, std::bind2nd(std::equal_to<char>(), '\t'));
		prot_id_clusters.push_back(vector<unsigned>());
		BOOST_FOREACH(const string& s, fields) 
		{
			unsigned prot_id_index;
			if (prot_id_strs.is_indexed(s, prot_id_index))
			{
				prot_id_clusters.back().push_back(prot_id_index);
				// cluster id is 0-based index 
				// (cluster 0 is empty and will be filled with unclustered items)
				prot_id_to_cluster_id[prot_id_index] =  prot_id_clusters.size() - 1;
			}
		}
	}
	args.VERBOSE(2) << "\t\t" << prot_id_clusters.size() << " predefined clusters\n";
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	debug_func_cerr(6);
	t_portable_timer timer;

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// saving parameters
	string exe_parameters;
	for (int i = 0; i < argc; ++i)
		exe_parameters += string(argv[i]) + ' ';

	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	// get command line arguments
	t_program_args args;
	try
	{
		if (!process_command_line_options(argc, argv, args))
			return 1;
	
		//
		//	parse pairs
		//
		//
		t_pairs				prot_pairs;
		t_index_string		prot_id_strs;
		do_parse_pairs parse_pairs(prot_id_strs,  prot_pairs, args.comment_letters);
		args.VERBOSE() << "\tReading in paired data...\n";
		seq_pairs::read_data(parse_pairs, args.file_input, args.VERBOSE());
		if (!prot_pairs.size())
			throw std::runtime_error("No pairs distances data supplied from STDIN.");
		args.VERBOSE(2) << "\t\t" << prot_pairs.size() << " pairs read.\n";
	






		
	
		//
		//	Read predefined clusters if necessary
		//
	
		// both predefined and single linkage clustered (but not parachuted)
		vector<unsigned>	prot_id_to_cluster_id(prot_id_strs.size(), 0);
	
		//
		// Which prot_id are in which cluster
		//
		// 		cluster 0 is for unclustered prot_ids
		//
		deque_vec_u			prot_id_clusters;
		prot_id_clusters.push_back(vector<unsigned>());
	
		if (args.file_cluster.length())
			read_predefined_clusters(prot_id_strs, prot_id_clusters, prot_id_to_cluster_id, args);
		
		t_hashuu		prot_id_to_parachute_id;
		deque_vec_u		parachute_id_to_cluster_ids;
		vector<char>	is_parachuted(prot_id_strs.size(), 0);
	
		
		//
		//	decide which cluster each pair goes into
		// 
		// 		"parachuted" prot ids don't start off with a cluster but end up linked with
		// 			all the clusters they are paired with
		//			& can therefore belong to > 1 cluster
		//
		if (prot_id_clusters.size() > 1)
		{
	
			// cluster predefined and parachute those which link into predefined
			cluster_by_predefined_clusters(	prot_pairs,
											prot_id_strs,
											prot_id_to_cluster_id,
											is_parachuted,
											prot_id_to_parachute_id,
											parachute_id_to_cluster_ids,
											prot_id_clusters,
											args);
	
#ifndef NDEBUG
//			cerr << string(50, '*') << "\n";
//			cerr << "After cluster by predefined (including parachuting)\n";
//			dump_id_pairs(prot_pairs, prot_id_strs);
#endif
			resolve_parachuted_pairs(		prot_pairs,
											parachute_id_to_cluster_ids,
											is_parachuted,
											prot_id_to_parachute_id,
											args);

#ifndef NDEBUG
//			cerr << string(50, '*') << "\n";
//			cerr << "After resolving parachuted pairs\n";
//			dump_id_pairs(prot_pairs, prot_id_strs);
#endif

	
	
			// sort by cluster
			sort(prot_pairs.begin(), prot_pairs.end());
		}
	
	
		// 
		// Single linkage cluster all pairs not in predefined clusters
		// 
		//  
		single_linkage_cluster(prot_pairs, prot_id_clusters, prot_id_to_cluster_id, args);
	
		
		args.VERBOSE(1) << "\t" << (prot_id_clusters.size() - 1) << " clusters in total\n";
	
		
		// write to file
		ofstream strm_summary;
		if (!args.results_summary.empty())
		{
			strm_summary.open(args.results_summary.c_str());
			if (!strm_summary)
				throw runtime_error("Could not open " + args.results_summary);
		}

		//
		// get sizes
		//
		vector<unsigned>	cluster_id_to_cluster_size(prot_id_clusters.size(), 0);
		unsigned cnt_output_clusters = 0;
		// first cluster is for unclustered pairs. Ignore
		for (unsigned i = 1; i < prot_id_clusters.size(); ++i)
			if (prot_id_clusters[i].size() >= args.minimum_cluster_size)
			{
				cluster_id_to_cluster_size[i] = prot_id_clusters[i].size();
				++cnt_output_clusters;
			}
		// get largest 10 sizes
		vector<unsigned>	largest_cluster_sizes(10);
		largest_cluster_sizes.erase(partial_sort_copy(	cluster_id_to_cluster_size.begin(), 
														cluster_id_to_cluster_size.end(),
														largest_cluster_sizes.begin(), 
														largest_cluster_sizes.end(),
														std::greater<unsigned>()),
														largest_cluster_sizes.end());

		largest_cluster_sizes.erase(std::remove(largest_cluster_sizes.begin(), 
												largest_cluster_sizes.end(), 0U), 
												largest_cluster_sizes.end()); 
		
		unsigned cnt_above_threshold_identifiers = get_identifier_count(	prot_id_strs,
																prot_id_clusters,
																cluster_id_to_cluster_size,
																args);


		// 
		// only print out identifiers one line per identifier
		// 
		//  
		if (args.identifiers_only)
		{
			sort_and_print_identifiers(prot_id_strs, prot_id_clusters, cluster_id_to_cluster_size, args);
			args.VERBOSE(2) << "\t" << string(60, '_') << "\n";
			string total_time = timer.pretty_reset();
			print_summary(	prot_id_strs, prot_pairs, prot_id_clusters, largest_cluster_sizes, 
							cnt_output_clusters, total_time, cnt_above_threshold_identifiers, 
							strm_summary);
			print_summary(	prot_id_strs, prot_pairs, prot_id_clusters, largest_cluster_sizes, 
							cnt_output_clusters, total_time, cnt_above_threshold_identifiers, 
							args.VERBOSE());
			return 0;
		}

		
		// 
		// annotate unclustered prot pairs with clusters
		// remember cluster 0 = orphan
		//  
		assign_unclustered_prot_pairs_to_clusters(prot_pairs, prot_id_to_cluster_id);


		// 
		// prot pairs belonging to small clusters can be marked for removal
		// remember cluster max_size should be removed
		//  
		args.VERBOSE(2) << "\tSorting " << prot_pairs.size() << " pairs by clusters \n";
		sort(prot_pairs.begin(), prot_pairs.end());

		if (args.minimum_cluster_size)
		{
			args.VERBOSE(2) << "\tMarking pairs in small clusters for removal\n";
			mark_prot_pairs_of_small_clusters_for_removal(prot_pairs, cluster_id_to_cluster_size);
			removed_marked_pairs(prot_pairs, args);
		}
		


		//
		print_pairs(prot_pairs,
					prot_id_strs,
					is_parachuted,
					prot_id_to_parachute_id,
					parachute_id_to_cluster_ids,
					cluster_id_to_cluster_size,
					args);
		args.VERBOSE(2) << "\t" << string(60, '_') << "\n";
		string total_time = timer.pretty_reset();
		print_summary(prot_id_strs, prot_pairs, prot_id_clusters, largest_cluster_sizes, 
					  cnt_output_clusters, total_time, cnt_above_threshold_identifiers, args.VERBOSE());
		print_summary(prot_id_strs, prot_pairs, prot_id_clusters, largest_cluster_sizes, 
					  cnt_output_clusters, total_time, cnt_above_threshold_identifiers, strm_summary);
	
	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}


